A simple DEB-based ecosystem model

Abstract A minimum stoichiometric carbon and nitrogen model of an entire ecosystem based on Dynamic Energy Budget (DEB) theory is presented. The ecosystem contains nutrients, producers, consumers, decomposers and detritus. All three living groups consist of somatic structure and either one (consumers and decomposers) or two (producers) reserve compartments, hence the living matter is described by seven state variables. Four types of detritus are distinguished. As the system is closed for matter, the dynamics of the nutrients carbon dioxide and ammonium follow automatically from the dynamics of the other 11 state variables. All DEB organisms in the model are V1-morphs, which means that surface area of each organism is proportional to volume. The resulting ontogenetic symmetry implies that complicated modelling of size structure is not required. The DEB V1-morph model is explained in detail, and the same holds for the idea of synthesizing units, which plays a key role in DEB modelling. First results of system dynamics are presented.


Introduction
Human activities often impose direct effects on organisms of high conservation value. The construction of wind farms, for example, puts iconic bird species, such as vultures, eagles and larger seabirds, at risks of collisions with wind turbines. Mortality rates increase, with possible negative effects on population size (Schippers et al., 2020). More subtle indirect effects through changes in the entire food web should, however, not be underestimated. Offshore wind farms in the southern North Sea decrease the sea surface wind speed and affect the stratification development (Christiansen et al., 2022). These changes in the physical environment may alter primary production level with cascading effects through the entire food web, ultimately up to seabird abundance. The use of mathematical ecosystem models is required for the assessment of such indirect effects. This, of course, holds for many other types of human impact, such as eutrophication, pollution or climate change-related seawater warming and acidification.
2008) provides models of organism's metabolism, which are now widely used at the level of the individual. The standard DEB model, appropriate for species with isomorphic growth, has been tested for over a thousand different species (Marques et al., 2018, van der Meer et al., 2014. The V1-morph DEB model, which is suitable for organisms growing in such way that surface area is more or less proportional to mass, has been used for a variety of micro-organisms (Brandt et al., 2004, Grossowicz et al., 2017, Lika & Papadakis, 2009, Livanou et al., 2019. Despite the availability of such a well-established theoretical framework, most community and ecosystem models still use descriptions of species metabolism that are not derived from a general theory, but are merely based on old habits (Anderson & Mitra, 2010, Babel et al., 2019, Maps & Record, 2020. Such practice hampers generality and transferability and thus slows down scientific progress. In this paper we describe in detail the simplest possible ecosystem (or community, if you like) model that is based on DEB theory. The model should be considered as a basic model that can be extended for specific conservation goals. This minimum ecosystem model includes a nutrient, a producer, a consumer and a decomposer, and has a closed mass-balance, essential prerequisites of ecosystem models. It is assumed that all species are V1-morphs, which means that area is proportional to volume for an individual organism. DEB theory assumes that some process rates, such as food ingestion rate, are area related, whereas other rates, like maintenance rate, are volume related. But for a V1-morph the areavolume ratio does not change during growth and all rates are therefore implicitly volume-related. Organism size plays no role, a condition sometimes referred to as ontogenetic symmetry (Roos et al., 2013). Take, for example, an organism in the form of a rope that only grows in length. It then does not matter whether the rope is cut in smaller pieces or not, as long the area of the tips is negligibly small. Hence, there is no real difference between individuals and populations and complicated modelling of the size structure of populations is not needed.
We embroider on the work of Kooijman & Nisbet (2000) and Kooijman (2010), who presented a first version of this model, which they called the canonical community model, although they did this in a rather concise way. We extended the model by including the possibility that the producer can also shrink when nutrient or light levels get too low. The main value of our contribution is that we provide an extensive treatment of the underlying ideas and how they translate in a set of differential equations. We also pay attention to parameter estimation and implementation.
Before we discuss the ecosystem model in detail, we first introduce the specifics of the V1-morph model and its link with the standard DEB model for isomorphs, discuss the use of different dimension frameworks (e.g. energy-length or mass-mass) in DEB modelling and explain the idea of synthesizing units (SUs), which play a key role in DEB theory. We assume some basic knowledge of the so-called standard DEB model (Kooijman, 2010), for which gentle introductions are available (van der Meer, 2006b(van der Meer, , 2016(van der Meer, , 2019.

V1-morphs
One of the main assumptions in DEB theory is that ingestion and assimilation rates are proportional to the surface area of the structural body of the organism. The standard DEB model deals with isomorphic animals. An isomorphic animal does not change in shape during growth. If two animals have the same shape, then every kind of length measure taken at one individual relates to the same measure taken at the other individual by a constant factor. Isomorphism thus implies that surface area scales with squared length and volume with cubic length. Although most animal species are approximately isomorphic, many other organisms are not. Filamentous fungi form branches, where the uptake of resources occurs at the tip area of each branch and the number of branches is proportional to the volume of the total mycelium (Gow & Gadd, 1995, Trinci, 1974. Such an organism is called a V1morph in DEB theory, which implies that the scaling factor that links area to volume equals one. Area is thus proportional to volume. Other (theoretical) examples of V1-morphs are organisms in the form of a rope that only grow in length and where the tips do not contribute to the uptake, or organisms in the form of sheets or crusts that only grow in diameter and where the uptake is only at the top or bottom (Fig. 1). Colonies of archaea, cyanobacteria or green algae may also fit this pattern. One specific example is the genus Pediastrum.
Organisms that do not adhere to isomorphs can be modelled using a shape correction factor. The shape correction factor is defined as where V is structural volume, V d is a reference volume required to keep the correct physical dimensions and where the scaling factor x indicates the morph type, having a value of 1 for a V1-morph, 2/3 for an isomorph and 0 for a V0morph, the latter thus having a constant surface area. The assimilation rate of a V1-morph equalṡ which can be simplified tȯ The parameter {ṗ Am } is one of the primary parameters of the standard DEB model and stands for the maximum areaspecific assimilation rate; the parameter [ṗ Am ] is the maximum volume-specific assimilation rate.

V1-morph reserve dynamics and structural growth
One consequence of the assumptions of homeostasis that are made in DEB theory ( for a detailed explanation, see Kooijman, 2010) is that reserve density [E], which is the amount of reserves per unit of structural volume, follows first order dynamics. In the case of V1-morphs, the rate at which the reserve density drops down when no assimilation takes place only depends upon the reserve density whereṗ A is the assimilation rate and f is the scaled functional response indicating the food availability. The proportionality coefficientk E has been given the name 'specific energy conductance' and has the physical dimension 'per time'. The maximum reserve density equals the ratio of the maximum volume-specific assimilation rate to the specific energy con- The reserves E equal the reserve density [E] times the volume V. Using equation 1 and the chain rule gives the mobilization rate, which is the rate at which the reserves are mobilized:ṗ The allocation of the mobilized matter can be simplified if organisms do not produce gonads, but simply divide into two daughter cells. Such dividing organisms are classified as juveniles in DEB terminology. For juveniles, maturity and maturity maintenance are not directly relevant for describing growth. This is a consequence of the κ-rule of DEB theory, which says that a fixed fraction κ of the mobilization rate is spent on growth and somatic maintenance. The remaining fraction is spent on maturity and maturity maintenance. To keep things as simple as possible, we therefore take κ = 1. Surface area-related maintenance costs are irrelevant for V1-morphs, and hence the allocation is given bẏ where [E G ] are the energetic growth costs per unit of growth of structural volume and [ṗ M ] are the volume-specific maintenance costs. Substituting equation 2 in equation 3 gives the growth equation Under constant food conditions, when the reserve density is in equilibrium and proportional to the scaled functional response, the growth equation simplifies to whereṙ is the specific growth rate. Hence, the growth rate of the structural volume is proportional to the structural volume itself. V1-morphs show exponential growth at constant food density.
The rate at which the reserves are mobilized is thuṡ

A mass-mass framework
So far the V1-morph DEB model has been written in a socalled energy-length framework, but, following Kooijman & Nisbet (2000), the model can easily be re-written in a massmass framework (Table 1).
The structural body can be expressed in terms of its mass, or better in terms of the amount of C-atoms, by is the specific density of the structural body expressed in C-moles per volume. In accordance with the strong homeostasis assumption [M V ] is a constant. Reserve density can be expressed as the number of C-atoms in the reserve per C-atom in the structural body. Reserve mass is related to reserve energy by M E = μ −1 E E, where μ E is the potential energy of the reserves expressed in energy per C-mole. Reserve mass relative to structural mass thus  Table 1: Primary parameters of the DEB model for V1-morphs, specific for a mass-mass framework. The last column indicates the relationship with the parameters from the energy-length framework that has been replaced.

Symbol
Dimension Interpretation Relationship

Fig. 2:
The preference SU. The unit will dissociate a less-preferred object (blue diamond) at the arrival of a preferred object (red diamond), but no product (yellow diamond) is formed at this transition from S2 to S1.
Within a massmass framework the DEB equation for reserve dynamics of V1-morphs thus looks like which can be re-written as where j EAm is the maximum mass-specific assimilation rate expressed in C-mole reserve per C-mole structure per time. Growth is given by where −1 refers to the mass-specific maintenance costs, expressed in C-mole reserve per C-mole structure per time, and gives the mass-specific costs of growth, expressed in C-mole reserve per C-mole structure. The growth equation has the form dM V /dt =ṙM V , whereṙ is the specific growth rate. Reserves are mobilized at a rate (k E −ṙ)m E M V .

Fig. 3:
A demand-driven preference SU has a constant product formation rate that equals the demand rate (black horizontal line), but the rate of product formation on the basis of the preferred substrate (red line) depends on the arrival rate of the preferred substrateJ 1 , horizontal axis. The difference between the black and red lines gives the rate of product formation on the basis of the less-preferred substrate, whose arrival rate is constant and equal as or larger than the demand rate. The blue line gives the rejected flux of preferred substrate. At high arrival rate of the preferred substrate hardly any less-preferred substrate is used.

Synthesizing units
DEB theory makes use of so-called synthesizing units or SUs. The idea of SUs could be formalized as a continuous-time Markov process, which is a stochastic process having the property that the future state of the unit (for example empty or busy) only depends upon the present and is independent of the past (Ross, 1989). Objects arrive at an SU by means of a Poisson process, which means that the number of arrivals in any interval is proportional to the length of the interval, which is equivalent to the assumption that the interarrival time is exponentially distributed (Ross, 1989). When the unit is empty, an arriving object is accepted and used by the unit for the synthesis of some product. But arriving objects are rejected in case the unit is busy with handling a previously arrived object. The time after arrival that it takes the unit to produce the product can be described as a stochastic process too, but might also be defined as a fixed time period. Different types of objects can be involved and these types are either substitutable, which means that none of them is essential for product formation, or non-substitutable, which means that all types are needed in the production process. The ecosystem model presented here uses several types of SUs. The first SU type handles one or more substitutable objects, which are some type of substrate. The second type receives two substrates, which are also substitutable, but in this case the unit has a preference for one substrate above the other. A specific type of this preference unit concerns the case when the rate of product formation is demand-driven, which means that it is constrained at some constant value. The third SU type receives two or more non-substitutable substrates, which implies that all types are required for product formation.

SU with one or more substitutable substrates
First consider an SU that only uses one type of substrate. The unit is either empty (a fraction S 0 /S of all units is empty), or contains a substrate (fraction S 1 /S), where S is the overall density of units. If a unit is empty it will accept an arriving substrate. The following two differential equations describe the dynamics of the number of units in each of the two states: whereJ i =ḃX is the arrival rate of substrates andk is the dissociation parameter. These differential equations are equal to zero at equilibrium, sometimes called pseudo-equilibrium because X is assumed constant. This reveals S * 1 =ḃ k XS * 0 , where S * 0 , and S * 1 refer to the equilibrium densities. Since The overall processing rate per SU equals the dissociation parameter times the fraction of units of type 1, and is thus given bẏ Ecologists will recognize that the result is equivalent to Holling's disc equation, also called the type II functional response equation, whereḃ is the area of discovery, and 1/k the handling time. The equation is often written aṡ wherek is the maximum processing rate and where the half- Similarly, an SU that uses multiple types of substitutable substrates results in the equivalence of Holling's type II functional response equation for multiple prey. For two substrates, the unit is either empty, contains a substrate of type 1 or of type 2. The three differential equations that describe the dynamics of the density of units in each of the three states are whereJ i =ḃ i X i is the arrival rate andk the dissociation parameter for substrate i. At equilibrium and since S * 0 +S * 1 +S * 2 = S, it follows that the overall processing rate for substrate i is given bẏ

Preference SU with two substitutable substrates
Consider an SU that can use two types of substrates, but has a preference for one type. The unit is either empty (a fraction S 0 /S of all units is empty), contains the preferred substrate (fraction S 1 /S) or the less appreciated substrate (fraction S 2 /S). If the unit is empty it will always accept both types once they arrive. If the unit contains the less appreciated substrate 2 and a preferred substrate 1 arrives, it will delete type 2 and replace it by the preferred type 1 (Fig.2) following differential equations describe the system: whereJ i =ḃ i X i is the arrival rate andk i is the dissociation parameter for substrate i. The dynamics of the preferred substrate 1 follow exactly Holling's disc equation, as it does not matter whether the unit is in state 0 or state 2, in both cases it will accept an arriving item of type 1. In equilibrium Noting that S 0 = S − S 1 − S 2 and setting the third differential equation equal to zero, giveṡ from which the equilibrium density for state 2 follows For each unit, the product P is delivered with a ratė where y i gives the number of product particles formed per substrate particle i = 1, 2.
In case of 'demand kinetics' the product delivery rateJ P has to be constantJ P =k P and the dissociation parametersk 1 anḋ k 2 will be tuned accordingly. The only requirement is that the ratio between the two parameters remains constant, given by the so-called preference parameter ρ =k 2 /k 1 (Fig.3).
In equilibrium, usingk 1 =k andk 2 = ρk, the delivery rate is given by and the dissociation parameterk can be solved from this equation, which can be re-written as a quadratic equation of and where the solution, which always has to be positive, equalsk An equivalent approach is to assume thatk 1 =k P /θ and to solve the equation (2010) uses still another approach, with x = S * 1 /(y 2 ρS * 2 ) as the unknown of a quadratic equation. All three approaches reveal the same result. Below, the first approach will be followed.

SU with two or more non-substitutable substrates
Consider, for example, an SU that requires two substrates, one of type A and one of type B, to produce a product P. Suppose further that the binding of one type does not interfere with that of the other, which may be called parallel processing. The density of empty units is given by S 00 , the density of those with only A bound by S 10 , with only B by S 01 and if both A en B are bound by S 11 (Fig. 4) In parallel processing of two non-substitutable substrates, the binding of one type of substrate does not interfere with that of the other. An SU that requires two substrates binds these in a random order, depending upon the order of arrival. When both are bound, the busy unit 11 will produce a product and then return to the empty state 00. apply: In equilibrium, As the equilibrium densities for all states of the unit are expressed in terms of the density of the busy unit, it is easy to arrive at the process rate per unit (equivalent to the production rate of product C), which, as before, equalsk times the fraction of busy unitṡ This expression can (after some algebraic manipulation) be simplified toJ whereJ Cm =k is the maximum process rate andJ A =ḃ A X A andJ B =ḃ B X A are the rates at which substrates A and B, respectively, arrive at each server. Thus,J −1 Cm is the expected 'servicing' time andJ −1 A andJ −1 B are the expected interarrival times of type A and B substrates, respectively.
It can be shown that if both types of substrates can not be bound simultaneously (i.e. the binding of one can only start when the other is already bound) and when the order in which they arrive is important, say first A and then B, (this is called sequential processing) that the process rate decreases tȯ The expected processing timeJ −1 C is then simply the sum of the expected 'servicing' timeJ −1 Cm and the expected interarrival times of the substrates (hereJ −1 The parallel processing SU described above behaves very much like a minimum operator (Fig. 5), where it is assumed, following Liebig's law, that the processing rate is only limited by one type of substrate, the so-called limiting substratė In ecology, the use of the minimum operator in growth models has been popularized by Tilman (1988  Contours of the production rate as a function of the rates at which substrate A and B arrive, for a parallel processing SU that requires two non-substitutable substrates to produce a product C (as in Fig. 4).
greater realism, but also because it prevents numerical problems related to the stepwise change of the minimum operator, in further applications, for example in population models.

The canonical ecosystem model
The canonical ecosystem model described here is the most simple DEB-based model for V1-morphs, formulated in a mass-mass framework and keeping track of carbon and nitrogen fluxes. It contains apart from the nutrients carbon dioxide and ammonia, three living (producers, consumers and decomposers) and one non-living group (detritus). Each living group is characterized by a structure and by two (producers) or one (the other two groups) type of reserves. Detritus is split up in four types, depending on the origin, which can either be consumer faeces as a result of eating producers, consumer faeces from eating decomposers, consumer dead structure or consumer dead reserves. The model has therefore 13 state variables ( Table 2). As the system is closed for matter, the dynamics of the nutrients follow automatically from the dynamics of the other 11 state variables. Fluxes from and to these 11 state variables are presented in a flux matrix (Table  3). Parameters and conversion coefficients are not explained in the text, but given in Tables 4, 5 and 6.

Consumers
Consumers predate on both producers and decomposers and the overall predation rate is given by the processing rate of an SU with two substitutable substrates, or in the terminology of an animal ecologist by a Holling's type II functional response with two prey types. The consumption or removal rate of producer structure (in terms of a flux measured in C-moles) is given byJ where the minus sign indicates that matter is subtracted from a compartment, in this case from producer structure.
All fluxes in the model that are subtracted from one of the 11 compartments have a negative sign, and those that are directed to one of these compartments have a positive sign. The removal of producers by predating consumers or the consumption of detritus by decomposers, for example, are represented by negative fluxes. The same holds for fluxes that describe the mobilization of reserves or structure. Fluxes directed towards reserves (assimilation) or towards structure (growth) have a positive sign. But beware that a general description of a flux, for example one indicating the required maintenance costs, without being specific where the flux comes from or goes to, is given as a positive figure.
The consumption of producer structure goes hand-in-hand with the consumption of the two producer reserve typeṡ for i ∈ {1, 2}. Part of the consumed matter ends up in the reserves of the consumeṙ Note that the minus sign will turn the negative consumption fluxes into a positive assimilation flux. Another part of the consumed matter ends up in the producer faeces fraction of the detritus poolJ Similarly for the consumption of decomposers, which only have one type of reserves, by consumerṡ The consumers mobilize their reserves at a rate (negative sign) equal toJ where the termṙ VC,GC m EC M VC is needed to avoid dilution by growth of the reserve density (compare with equation 6).
The instantaneous growth rate equalṡ where j EC,MC is the mass-specific maintenance rate of the consumers (see also equation 8).
Usually, the mobilized energy from the reserves is first spent on maintenance, and the rest goes to growth. The maintenance costs, when paid from the reserves, equalJ EC,MC = j EC,MC M VC . However, it might be that the mobilized flux from the reserves is not sufficient to pay the required maintenance. It is therefore assumed that maintenance can also be paid from a flux that comes from structure, and thus allowing shrinkage of structure. The maintenance costs, when paid from structure equalJ VC,MC = j VC,MC M VC . When reserves are transformed into structure, it is assumed that y EC,VC C-moles of reserve are needed to create one C-mole of structure. For convenience, it is also assumed that maintenance, when paid from reserves, requires a fraction y EC,VC more C-moles than when it is paid from structure. Hence, the mobilized flux from structure can be expressed in terms of C-moles of reserves by multiplication with y EC,VC .
It is assumed that the two fluxes, one containing reserve substrate and the other structure substrate, are directed towards a demand-driven preference SU that has a preference for reserves above structure. The demand is equivalent to the overall maintenance costs and equals, in terms of reserve C-moles,J EC,MC = j EC,MC M VC . The mobilized flux from structure can for convenience be taken just as large as to be able to pay for all maintenance. It is, in terms of C-moles of structure, given byJ It is just sufficient for paying all maintenance, which implies that The work of the preference SU results in a rate of maintenance costs paid from reserves (using the notation of the preference SU as used in Preference SU section and given as a positive flux) equal toJ    where the unknown dissociation parameter is given bẏ with the coefficients A, B and C as given earlier in Preference SU section, and withJ 1 = −J EC,CC ,J 2 = −y EC,VCJVC,CC , k p =J EC,MC , y 1 = 1 and y 2 = 1.
For completeness, the equations can be written out explicitly. The dissipation flux from the reserves, which is equivalent to minus the rate of maintenance costs paid from reserves, is thus given byJ Half-saturation constant of light assimilation by producer reserve 1 Half-saturation constant of light assimilation of producer reserve 2  Martiny et al. (2014) n N,ED mole N/mole C 0.15 Decomposer reserves Martiny et al. (2014) withk as above and where ρ is a preference parameter.
Note that the overall maintenance costs are the sum of the costs paid from reserves and those paid from structure, which then directly gives the maintenance costs paid from structure, in terms of C-moles reserve, Recall that all these maintenance costs,J MC VC ,J MC EC andJ EC,MC , are presented as positive fluxes. Consequently, the dissipation flux from structure, which is equivalent to the part paid from structure, equals, in C-moles of structurė Alternatively, the maintenance costs paid from structure could have been calculated (using the notation of the preference unit (k +J 1 )(ρk +J 1 +J 2 ) but this gives, of course, exactly the same result.
The rejected flux from the structure (which equals -J VC,CC −J MC VC ) is immediately channeled back to the structure. The rejected flux from the reserves is allocated to growth, and the flux from reserves to growth therefore equalṡ which results in an increase of structure equal tȯ The overall mass balance of the structure is thus given by a So in fact knowingJ EC,CC andJ EC,MC is sufficient for describing the mass balance of the structure. One might therefore argue that the idea of a preference SU is not really needed. One simply can accept that a negative growth flux from structure to reserves is possible. However, for producers, which are discussed below and which have multiple reserves, such simplification is not possible and preference SUs are needed to describe maintenance and growth.
As shown above, consumers will not die from starvation, for example when maintenance costs cannot be met as a result of too low reserve density, but will then also use structure to pay for maintenance and consequently shrink. In addition, we assumed a natural background mortality, where the death rate depends upon the reserve density. The death of consumers results in a flux from consumer structure to the dead consumer structure part of the detritus pool, given bẏ J PV,HC = −J VC,HC =ḣ a M VC m EC /y EC,VC 1 + m EC /y EC,VC and in a similar flux of consumer reserves to the dead consumer reserve part of the detritus pool J PE,HC = −J EC,HC = m ECJPV,HC .

Producers
The producers have two types of reserves. The first type contains only carbohydrates and the assimilation rate for this type depends upon the light conditions and the carbon-compound concentration. Using an SU with two non-substitutable substrates gives the assimilation ratė where The assimilation rate of the other reserve type depends upon light, the carboncompound and ammonium and is (using a SU with three nonsubstitutable substrates) given bẏ where From each of the two reserves, a mobilized flux, which equalsJ is, together with a mobilized flux from the structure, which expressed in C-moles reserve equals y E i P,VPJV i P,CP = −y E i P,VP j V i P,CP M VP directed towards a demand-driven preference SU. So there are two of these SUs. Each of these two units generates a maintenance flux with contributions from both reserves and structureJ

13
The maintenance that is paid from the reservesJ MP E i P is calculated in a similar way as for the consumers (i.e. as explained above). This also holds for the maintenance paid from structure, which is given by y E i P,VPJ MP V i P . The rejected fluxes from the mobilized fluxes from structure are directly channeled back to the structure. Hence, these fluxes do not play a role in the mass balance. The two rejected fluxes from the mobilized fluxes from the reserves are now first directed towards a growth SU, which is a SU at which two non-substitutable substrates arrive. These two fluxes equal, in terms of C-moles of structure, The product of this growth SU represents the overall growth flux and is given bẏ The overall growth rate equalṡ Note that the three terms in the numerator at the right side of the equation depend upon the growth rateṙ VP,GP as well and no explicit equation for the growth rateṙ VP,GP is available. A root-finding procedure has to be used to calculate the growth rate.
The flux towards growth thus equalṡ and the contribution of each of the two reserves iṡ Recall that the two types of reserves arriving at the growth SU are non-substitutable substrates, and the substrate that is non-limiting at arrival will be rejected. A fraction κ E i of the rejected flux is channeled back into the reserves, the rest is dissipated. Each rejected flux equalṡ SUs are presented by circles. Nutrients and light arrive at two assimilation SUs, whose products enter the reserves (product flows are shown by blue arrows). For each of the reserves, the mobilized flux is, together with a flux from the structure, directed towards a maintenance SU. The products of these SUs are dissipated. The rejected fluxes are shown by red arrows. Those originating from the structural volume are sent back, and those originating from the reserves go to a growth SU, whose product is channeled to the structural volume. Fluxes rejected by the growth unit are partly sent back to the reserves, the rest dissipates. Nutrient fluxes rejected by the assimilation SUs are not depicted.
The total dissipation flux from the reserves, which also includes the maintenance paid from the reserves, is therefore equal toJ The dissipation flux from the structure equals the maintenance costs paid from the structure (in units of structure) Summarizing, nutrients and light arrive at two SUs that handle non-substitutable substrates. The products are channeled to two reserve types. For each of the reserves, the mobilized flux is, together with a flux from the structure, directed towards a demand driven preference SU, whose product is used for maintenance. The rejected fluxes originating from the two reserves go to a growth SU, which handles two nonsubstitutable substrates. All together, five SUs are thus needed to describe the mass fluxes through the producer (Fig. 6).

Decomposers
The decomposers live on detritus and the removal rate for each of the four types of detritus follows a multiple prey type

Temperature and light
DEB theory assumes that all rates respond in the same way to temperature (Kooijman, 2000). The Van't Hoff-Arrhenius equation is used to correct rates for temperaturė wherek(T) is the reaction rate that depends upon the absolute temperature T (in Kelvin),k 1 is the reaction rate at a reference temperature T 1 and T A is the Arrhenius temperature, assumed to be constant within the same organism. Temperature corrections can be added to the main equations by multiplying all rates by the Van't Hoff-Arrhenius correction factor.
In the canonical ecosystem model, seasonal fluctuation in temperature is modelled with a cosine function that varies between the minimum temperature T min ( • C) and maximum temperature T max ( • C), which is reached at day of the year t Tmax : The constant 273.15 converts temperature in degree Celsius to Kelvin, as required by the Van 't Hoff-Arrhenius function.
Yearly fluctuation in light irradiation (J L in in MJ m −2 day −1 ) is modelled similar to yearly fluctuation in temperature, but using different parameters:

15
We assumed that light availability decreases through self shading as a function of producer structural mass. Light availability including self-shading is given bẏ Parameter K S modifies the strength of self-shading by the producer.

Parameterization
All model parameters and their sources are listed in Tables 5  and 6. Where possible we refined the original parameterization of Kooijman & Nisbet (2000). Producer parameters were mainly taken from Lorena et al. (2010), who proposed a DEB model for a general microalgae. Producer specific energy conductance parametersk E 1 P andk E 2 P were refined by Grossowicz et al. (2017). Parameters for the decomposer were taken from Eichinger et al. (2009), who presented a DEB model for marine bacterial biomass. Parameters related to the decomposer's functional response were taken from Lorkowski et al. (2012).
Consumers are for reasons of convenience treated as V1morphs in the canonical ecosystem model, although most zooplankton and zoobenthos species are actually better characterized as isomorphs. V1-morph parameters can be derived from isomorph parameters by assuming that all individuals have the same reference structural volume V d . The maximum volume-specific assimilation rate can then be obtained from the maximum area-specific assimilation rate as follows: Similarly, the specific energy conductance parameterk E can be obtained as ki E = viV d −1/3 , wherė v (cm d −1 ) is the primary DEB parameter that describes the energy conductance of an isomorph. Subsequently, the energylength parameters for V1-morphs can be converted to massmass parameters by dividing each volume-specific parameter by the energy density of structure μ E [M V ] in J cm −3 .
We derived V1 parameters for consumers from the isomorphic parameters for Daphnia magna available in the Add-my-Pet library (Kooijman & Gergs, 2019). We choose the length at puberty for D. magna as reference structural length (L d = V 1/3 d = 0.07 cm) and a value of μ E [M V ] equal to 3556 J cm −3 (Table 4). The obtained mass-specific maximum assimilation was divided by the estimated digestion efficiency of food to reserve (κ X = 0.9) to obtain a massspecific maximum consumption rate (j EAm ). We assumed that in the canonical ecosystem model, maximum consumption of producer structure equals j EAm , while maximum consumption of decomposer structure by consumers equals 0.5j EAm . Consumer mortality was set at h C = 0.01.
Parameters that describe the nitrogen content per carbon were taken as much as possible equal to the observed global median C:N ratio of 6.6 ('Redfield ratio') as described for phytoplankton (Martiny et al., 2014). Nitrogen content of the carbon and the nitrogen reserves of the producer were set to 0 and 0.8, respectively.
Temperature parameters were obtained from data on water temperatures in the North Sea from 2004 onwards at the measurement station of Schiermonnikoog Noord, The Netherlands (Deltares, 2021). Irradiance parameters were derived from measured irradiance from 2004 onwards at De Kooy, a measurement station of the Royal Netherlands Meteorological Institute (KNMI, 2021).

Implementation
The model was implemented in R software (R Core Team, 2021) using the R-package deSolve (Soetaert et al., 2010). The code to run the analysis is available via the Open Science Framework (https://osf.io/nfqgh/).

Dynamics
The dynamics of the canonical ecosystem model converge towards a stable attractor with two distinct producer (phytoplankton) peaks per year; one in April/May and a slightly lower and shorter peak in August (Fig. 7). Both producer peaks are followed by an increase in consumer biomass. Because the growth of producer biomass has already ceased before consumer biomass increases, the producer growth is limited by nutrient and light availability, which decrease due to assimilation and self-shading. However, consumers are responsible for the subsequent decrease in producer biomass, because producer background mortality is not accounted for. Because all detritus pools are associated with the consumer either through feeding (faeces: M PP and M PD ) or dead consumer biomass (M PV and M PE ), detritus biomass increases alongside the growth of consumers. Upon suppression of producer biomass, consumers die from starvation and background mortality and the detritus is cleared by decomposers.

Discussion
We present a minimum ecosystem model based on DEB theory, with an example mimicking a marine ecosystem. The model is described in a mass-mass framework, as opposed to the more conventional energy-length DEB framework (Kooijman, 2010), and explicitly considers nitrogen and carbon stoichiometry. All three functional groups, that is producers, consumers and decomposers, are treated as V1morphs and either contain one (consumers and decomposer) or two (producers) reserve compartments. Different types of SUs regulate substrate uptake, maintenance and growth. For the marine ecosystem application the producer represents the community of phytoplankton species, the consumer represents the zooplankton species and the decomposer represents the bacteria. Dynamics are driven both by environmental variables and species interactions. The set of parameters and simulations was presented for a part of the North Sea ecosystem, along the Dutch coastline. Under a seasonally fluctuating temperature and light regime, model dynamics show a spring and autumn peak of primary production (Daewel et al., 2014, Maar et al., 2018. A bimodal seasonal pattern of phytoplankton has been observed in coastal waters such as the North Sea (Cadéee, 1986, Colebrook, 1982, 1984, Desmit et al., 2020, Pitois et al., 2012, Reid et al., 1990 and is reproduced by several high-resolution ecosystem models of lower-trophiclevel processes (Los et al., 2008). Simulated consumer biomass also shows two peaks directly following the peaks in primary production. Despite large spatio-temporal variation, seasonal dynamics of zooplankton abundance in the North Sea typically show a unimodal pattern with high zooplankton biomass from spring until autumn (Colebrook, 1984, Fransz et al., 1991, Pitois et al., 2012, Van de Wolfshaar et al., 2021, although localized bimodal patterns in zooplankton have been observed in the North Sea (Colebrook, 1984). A potential reason for the difference between generally observed and modelled zooplankton seasonal dynamics is the strong coupling between producers and consumers in our model. Further explorations should reveal whether a different parameter setting can dampen the strong producerconsumer interaction, or that additional mechanisms, such as for example consumer interference or extra producer or consumer functional groups are required for a more sustained peak of consumer biomass. Simulated decomposer biomass shows two delayed peaks, but data to confirm whether such pattern occurs in the field are lacking.
Because the presented ecosystem model is based on a formal theory about metabolic organization (DEB), modelled processes and their parameters directly relate to processes at the individual level. This approach differs fundamentally from most marine ecosystem models of lower trophic levels (for an overview, see Maar et al., 2018), which often rely on ad hoc descriptions of population or ecosystem processes (Marques et al., 2014). In contrast, DEB provides a quantitative approach built on a collection of assumptions and stylized facts concerning individual-level energy acquisition, allocation and use and is well-tested across a wide range of taxa (Lorena et al., 2010, 2010, van der Meer et al., 2014. The use of DEB for ecosystem modelling also allows extending our minimum ecosystem model with additional food web components, as organismal DEB models and their parameters are readily available (Marques et al., 2018, van der Meer, 2006a, van der Meer et al., 2014. We provide an example on how to relate 'energy-length' DEB parameters for isomorphs to the 'mass-mass' parameters for V1-morphs as used in the present ecosystem model. If the V1-morph approach is deemed inappropriate, which is probably the case for higher trophic levels, a cohort-based approach could be adopted (van der Meer, 2016). However, this would considerably complicate model implementation and analysis. Future studies should shed more light on how system dynamics is affected by the assumption of V1-morphy for organisms that show a much larger size range during ontogeny than isomorphs that simply divide into two daughter cells. In addition, a spatial component can easily be incorporated through a coupled hydrodynamic model that simulates the transport of nutrients and pelagic ecosystem components across distinct spatial compartments. As such, an extended version of the model presented here would provide a powerful tool for studying a wide variety of questions related to marine conservation and management, such as ecosystem effects of human impacts (e.g. eutrophication, fishing), potential for mariculture and marine spatial planning. The use of the massmass framework allows to use stoichiometry of the different food web components to keep track of and maintain mass balance of nutrients such as nitrogen and phosphorous. In the example presented here, phosphorous is not included, but it can be added to address questions on for example phosphorous and nitrogen limitation in response to (terrestrial) management actions.
In its current form, the present model sacrifices realism to precision and generality (Levins, 1966). This contrasts most ecosystem models of lower trophic levels, which have high spatial and temporal resolution and consider multiple functional phytoplankton and zooplankton groups in order to maximize correspondence with data and by doing so try to achieve high predictive ability and thus realism. Due to their long computation times, these models rarely address parameter uncertainties nor study how the obtained results depend on model assumptions (Rose et al., 2010). This can more readily be achieved by simpler and more general models, for which a more complete sensitivity analysis is within reach. In addition, models like the present one can be used to explore how various ecosystem processes, such as the mode of trophodynamic control or the processes that limit phytoplankton and zooplankton growth, depend on various model parameters and assumptions.